%% Generate random GP functions (two-dimensional)
% choose covar function
covfunc = 2;
switch covfunc
    case 1; 
        k = @(x,y) 1*x'*y;                      % Linear
    case 2;
        k = @(x,y) exp(-100*(x-y)'*(x-y));      % squared exp
    case 3;
        k = @(x,y) exp(-1*sqrt((x-y)'*(x-y)));  % Ornstein-Uhlenbeck

end

% Choose points at which to sample
points = (0:0.05:1)';
[U, V] = meshgrid(points,points);
x = [U(:) V(:)]';
n = size(x,2);

% construct the covariance matrix
C = zeros(n,n);
for i = 1 : n
    for j = 1 : n
        C(i,j) = k(x(:,i),x(:,j));
    end 
end

% Sample from a Gaussian process at these points
u = randn(n,1);                 % sample u ~ N(0,I)
[A,S,B] = svd(C);               % factor C = ASB*
z = A*sqrt(S)*u;                % z = A S^.5 u ~ N(0,C) 

% Display function
figure; clf;
Z = reshape(z,sqrt(n),sqrt(n));
surf(U,V,Z);